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NUMERICAL BIFURCATION ANALYSIS OF CONFORMAL 
FORMULATIONS OF THE EINSTEIN CONSTRAINTS 

M. HOLST AND V. KUNGURTSEV 

Abstract. The Einstein constraint equations have been the subject of study for 
more than fifty years. The introduction of the confornial method in the 1970's 
as a parameterization of initial data for the Einstein equations led to increased 
interest in the development of a complete solution theory for the constraints, with 
the theory for constant mean curvature (CMC) spatial slices and closed manifolds 
^^ completely developed by 1995. The first general non-CMC existence result was 

establish by Hoist et al. in 2008, with extensions to rough data by Hoist et al. 
in 2009, and to vacuum spacetimes by Maxwell in 2009. The non-CMC theory 
remains mostly open; moreover, recent work of Maxwell on specific symmetry 
models sheds light on fundamental non-uniqueness problems with the confor- 
>0 mal method as a parameterization in non-CMC settings. In parallel with these 

mathematical developments, computational physicists have uncovered surprising 
I— I behavior in numerical solutions to the extended conformal thin sandwich for- 

•^r mulation of the Einstein constraints. In particular, numerical evidence suggests 

the existence of multiple solutions with a quadratic fold, and a recent analysis 

of a simplified model supports this conclusion. In this article, we examine this 

apparent bifurcation phenomena in a methodical way, using modern techniques 

C^ in bifurcation theory and in numerical homotopy methods. We first review the 

H evidence for the presence of bifurcation in the Hamiltonian constraint in the time- 

'— ^ symmetric case. We give a brief introduction to the mathematical framework for 

^_ analyzing bifurcation phenomena, and then develop the main ideas behind the 

►j^ construction of numerical homotopy, or path-following, methods in the analysis of 

^SJ bifurcation phenomena. We then apply the continuation software package AUTO 

\C to this problem, and verify the presence of the fold with homotopy-based numer- 

^^ ical methods. We discuss these results and their physical significance, which lead 

^^ to some interesting remaining questions to investigate further. 
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1. Introduction 

Einstein's gravitational field equations for relating the space curvature at a time 
slice to the stress-energy can be split into a set of evolution and constraint equa- 
tions. The four constraint equations, known as the (scalar) Hamiltonian constraint 
and the (3- vector) momentum constraint, constrain the induced spatial metric gij 
and extrinsic curvature Kij. The Einstein constraint equations have been the sub- 
ject of study for more than fifty years (cf. [2]). The introduction of the conformal 
method in the 1970's as a way of parameterizing initial data to the Einstein equa- 
tions led to increased interest in the development of a complete solution theory 
for the constraints, with the theory for constant mean curvature (CMC) spatial 
slices and closed manifolds completely developed by 1995. The CMC theory on 
closed manifolds is particular satisfying, with nearly all physically interesting cases 
exhibiting both existence and uniqueness of solutions. 

However, other than the near-CMC result of Isenberg and Moncrief in 1996 [10], 
the theory for non-CMC solutions remained completely open until the first far- 
from-CMC existence result was establish by Hoist et al. [: ] in 2008, with extensions 
to rough data by Hoist et al. in 2009 ["] , and to vacuum spacetimes by Maxwell in 
2009 [13]. However, the non-CMC theory remains mostly open, and what is known 
is much less satisfying than the CMC case; the new non-CMC results of Hoist et al. 
and Maxwell are based on new types of topological fixed point arguments, and while 
they establish existence, these arguments do not give uniqueness. Moreover, more 
recent work of Maxwell on specific symmetry models show in fact that uniqueness 
is lost, and also sheds light on some fundamental problems with the conformal 
method itself as a parameterization of initial data on manifolds that are not very 
close to CMC. 

Several decompositions of the equations have been formulated in addition to the 
conformal method, although the solution theory for only the conformal method is 
well-developed (cf. ['>]). In particular, in the extended conformal thin sandwich 
(XCTS) decomposition of the constraints, a conformal factor %p, the lapse iV, and 
the shift /3j, are solved for in five coupled elliptic equations with supplied back- 
ground data. In parallel with the mathematical developments in the theory for the 
conformal method, computational physicists have uncovered surprising behavior in 
numerical solutions to the XCTS formulation. In particular, numerical evidence 
suggests the existence of multiple solutions; Pfeiffer and York in [11] numerically 
construct two solutions for a specific choice of free data. 

In addition, there have been more theoretical results suggesting non-uniqueness. 
Under additional assumptions of conformal flatness and time-symmetry, the Hamil- 
tonian constraint becomes a decoupled equation for the conformal factor. Baum- 
garte, Murchadha and Pfeiffer ['>] solve this form of the Hamiltonian constraints 
using Sobolev functions, to find that an equation parameter has two admissible 
values consistent with the constraints. Walsh [ ] performs a Lyapunov-Schmidt 
analysis (a standard technique in bifurcation theory) of this form of the Hamiltonian 
constraint, and shows that at a critical point with two expected solution branches, 
the solution branches should take the form of a quadratic fold. In this article, we 
examine this apparent bifurcation phenomena in a methodical way, using modern 
techniques in bifurcation theory and in numerical homotopy methods. 

Outline of the paper. The remainder of the paper is structured as follows. In §2, 
give an overview of the constraint equations in the Einstein system and in particular 
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the XCTS formulation thereof. In §3, we give a brief introduction to the mathemat- 
ical framework for analyzing nonlinear operators in general. We build on this in §4 
where we expound on the theoretical framework underlying bifurcation analysis. In 
§5, we develop the main ideas behind the construction of numerical homotopy, or 
path-following, methods in the numerical treatment of bifurcation phenomena, and 
subsequently in §6 present the methodology of the numerical techniques we perform 
for this problem. In §7, we apply the continuation software package AUTO to the 
constraint problem, and verify the presence of the fold with homotopy-based nu- 
merical methods. We confirm the earlier results, as well as provide a framework for 
a more careful exploration of the solution theory for various parameterizations of 
the constraint equations. Analyzing the Hamiltonian constraint for time-symmetric 
conformally fiat initial data, as in [' '•, "], we demonstrate the existence and location 
of a critical point, evidence that the solution branch at the critical point forms a 
one-dimensional fold, and the form of the solution as continued past the critical 
point. 

2. CONFORMAL THIN SaNDWICH DECOMPOSITION 

In General Relativity it is common to look at the curvature of spatial hyper- 
surfaces taken at time-slices of space-time. The Einstein constraint equations are 
conditions on the induced spatial metric gij and the second fundamental form Kij 
for being a spatial slice. In addition, there are six evolution equations that govern 
how this geometric data evolves in a full spacetime. In solving the equations, a 
variety of decompositions have been proposed. In the XCTS formulation, proposed 
by York [ij], the shift vector /3j, the lapse N and a conformal factor if) are solved 
for given a set of supplied data that includes [gij, Uij, K, dfK), where gij is the con- 
formally related induced spatial metric, Uij its time derivative, and K the trace of 
the extrinsic curvature. The lapse and shift are formed from taking a level set of t 
and looking at the normal to the hypersurface n = N~^{dt — /3*(9j) where i indexes 
over spatial coordinates [i]. The Hamiltonian constraint of the XCTS equations 
can be written as 

V^V' - Im - ^K^^' + l^-'A'^A,, + 27r^V = (2.1) 

o 12 o 

Here V represents the covariant derivative, R the trace of the Ricci tensor and A 
the trace free part of the extrinsic curvature, all associated with the conformally 
scaled metric cjij. Following [5, 16] we assume time-symmetry (so Kij = 0) and 
conformally fiat initial data. This constraint then reduces to 

V^V' + 27rpV^^ = (2.2) 

Following [5] we let p be the constant mass-density of a star. Without loss of 
generality, we take p = outside r = 1 and look for solutions to (2.2) with boundary 
conditions 

f =0. r^O (2.3) 

ip = l, r = \ (2.4) 

We can note that the maximum principle does not apply for this equation and 
hence uniqueness is not guaranteed. We will see that this equation has two, one, 
or no solutions, depending on the value of p. 



4 m. holst and v. kungurtsev 

3. Nonlinear Operators on Banach Spaces 

Let X and Y be Banach spaces, and let X' and Y' be their respective dual spaces. 
Given a (generally nonlinear) map F: X -^ Y, we are interested in the following 
general problem: 

Find ueX such that F{u) = OeY. (3.1) 

We will give a brief overview of techniques for analyzing solutions to (3.1), and for 
characterizing their behavior with respect to parameters. 

We say that the problem F{u) = is well-posed if there is (a) existence, (b) 
uniqueness, and (c) continuous dependence of the solution on the data of the prob- 
lem. Recall that if F is both one-to-one (injective) and onto {surjective), it is 
called a bijection, in which case the inverse mapping F~^ exists, and we would have 
both existence and uniqueness of solutions to the problem F{u) = 0. Recall that 
F: X ^ Y is a continuous map from the normed space X to the normed space 
Y if \im.j^ooUj = u implies that \im.j^Qo F{uj) = F{u), where {uj} is a sequence, 
Uj G X. If both F and F~^ are continuous, then F is called a homeomorphism. 
If both F and F~^ are differentiable (see the below for the definition of differenti- 
ation of abstract operators in Banach spaces), then F is called a diffeomorphism. 
If both F and F~^ are fc-times continuously differentiable, then F is called a C''- 
diffeomorphism. A linear map between two vector spaces is a type of homomorphism 
(structure-preserving map); a linear bijection is called an isomorphism. 

We will need to assemble just a few basic concepts involving general nonlinear 
maps in Banach spaces, to provide the mathematical framework for the discussions 
in the next section. In particular, the following notion of differentiation of maps on 
Banach spaces will be required (see [15, 17] for more complete discussions). 

Definition 3.1. Let X and Y he Banach spaces, let F : X ^ Y , and let D d X 
he an open set. Then the map F is called Frechet- or F-differentiable at u E D if 
there exists a hounded linear operator Fu{u): X ^Y such that: 

lim -——\\F{u + h)- F{x) - FMWWy = 0. 
11^11x^0 \\h\\x 

The bounded linear operator Fu{u) is called the F-derivative of F at u. The 
F-derivative of F at m can again be shown to be unique. If the F-derivative of F 
exists at all points u E D, then we say that F is F-differentiable on D. If in fact 
D = X, then we simply say that F is F-differentiable, and the derivative Fu{-) 
defines a map from X into the space of bounded linear maps, F„: X — )■ C{X,Y). 
In this case, we say that F G C^{X; Y). Many of the properties of the derivative of 
smooth functions over domains in M" carry over to this abstract setting, including 
the chain rule: If X, Y, and Z are Banach spaces, and if the maps F: X —> Y and 
G: Y -^ Z are differentiable, then the derivative of the composition map H = GoF 
also exists, and takes the form 

HM = {G o FUu) = Gf{F{u)) o F„(n), 

where H^: X ^ C{X,Z), F^: X ^ C{X,Y), and Gp: Y -> C{Y,Z). Higher 
order Frechet (and Gateaux) derivatives can be defined in the obvious way, giving 
rise to multilinear maps and giving meaning to the notation F G C'^{X; Y). Note 
that below we will often encounter functions of two variables F{u, A), and will be 
interested in the derivatives of such functions with respect to each variable; we will 
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denote these using the consistent notation F„ and Fx- See [1, 15] for more complete 
discussions of general maps and differentiation in Banach spaces. 

A fundamental concept concerning linear operators that we will need for the 
discussions below is that of a Fredholm operator. Let X and Y be Banach spaces, 
and let A e C{X, Y), or in other words, A is a bounded linear operator from X to Y. 
In the case that X and Y have additional Hilbert space structure, wherein there 
is an inner-product, the Riesz Representation Theorem implies that the adjoint 
operator A* G C(Y, X), defined as the operator for which {Ax, y) = (x, A*y), exists 
uniquely. There are four fundamental subspaces of X and Y associated with A and 
A*, namely: 

(i) J^{A) : null space (or kernel) of A 

(ii) TZiA) : range space of A 
(iii) Af{A*) : null space (or kernel) of A* 
(iv) TZ{A*) : range space of A* 

One can consider the dimension (dim) and co-dimension (codim) of each of these 
four spaces, where co-dimension is taken to be relative to the larger spaces that 
they are subspaces of. The operator A G C{X, Y) is called a Fredholm operator if 
and only if: 

(i) dim{Af{A)) < cx) 
(ii) dim{n{A)) < cx) 

The difference of these two dimensions, namely 

md{A) = dim{Af{A)) - dim(7^(A)) (3.2) 

is called the Fredholm index of A. A basic result about Fredholm operators is the 
following. 

Theorem 3.2. If A E C{X,Y) is a Fredholm operator, then 

(i) TZ{A) is closed. 

(ii) A* is Fredholm with index ind(y4*) = — ind(y4). 
(iii) K is a compact operator, then A+K is Fredholm with ind{A+K) = md{A). 

Proof. See [17]. D 

4. Bifurcation Theory for Nonlinear Operators Equations 

Bifurcation Theory studies the branching of solutions as governed by param- 
eter(s). Throughout this section we will refer to a "solution curve" which, for 
F{u,X) plots points in (||m||,A) space at which (m. A) solves F{u,X) = 0, and ||m|| 
denotes the norm of u. Specifically, we are interested in the local behavior of the 
solution curve in a neighborhood of a known solution (mq, Aq). This is because, in 
order to explore the solution space to a certain operator equation F{u, A), we often 
solve F{u, Aq) for u, then obtain another {ui, Ai) from the original solution data we 
obtain, and similarly continue to subsequent solutions. 

Bifurcation theory has its foundation in the implicit function theorem. 

Theorem 4.1 (Implicit Function Theorem). // it holds that, for an operator F : 
X X M. -^ Y , if F{u, A) with F{uq, Aq) = 0, F and F^ are continuous on some 
region U xV with (xo,Ao) G f/ x V^, and Fu{uo,Xo) is nonsingular with a bounded 
inverse, then there is a unique branch of solutions (m(A),A)) in for X E V, e.g. 
F[u{X), A) = 0. Furthermore u{X) is continuous with respect to X in V 
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This theorem states that if the operator F^ is nonsingular at a certain point 
(xo, Ao) there is a unique solution u for each A close to Aq on either side of Aq and 
we can plot a one-dimensional curve in (||m||,A) space through (||mo||,Ao). With a 
singular F„, however, such a branch is not guaranteed, suggesting the possibility 
of two or more such u{X) branches or no solutions for some A in every neighbor- 
hood around Aq- The exact form of the branching depends on the dimensions of 
Fu{xo, Ao), Fx{xo, Ao) and whether Fa(xo, Ao) G 7^(F„(xo, Aq)). 

In the case of a "fold", wherein there is still a one-dimensional path through 
(mq, Aq) but the path exists solely on one side of Aq for A, we have 

(i) dim(Ar(F,(Mo,Ao))) = l 
(ii) dim(A/'(-Pu('Uo5 Ao)*)) = 1, where F„(mo,Ao)* is the adjoint operator of 

Fu{Uq,\q) 

(iii) FA(Mo,Ao)^7^(F„(Mo,Ao). 

It can be shown, under these circumstances, near (mq, Aq) the continuation of the 
solution will be of the form 

u{e) = Mo + 60 + Cit^ + 0{t^) (4. 1) 

A(e) = Ao + C2t^ + 0(e3) (4.2) 

where is a basis for the null-space of Fu{uq,Xq) (see [ i]). Note that there is no 
linear A(e) term, so the solution curve, at first approximation, stays at a constant 
value of A and changes in u. The specific values of Ci and C2 are based on the 
Hessian of the operator F. With a fold the C2 will be negative, so the solution 
curve shows A increase up to a certain critical point Ac then decreases. With Ci < 
the actual solution u changes at first-order along 0, the change then dampens in 
magnitude. Hence with C2 7^ 0, the resulting solution curve appears locally as a 
sideways parabola, and hence called a "simple quadratic fold". If, on the other 
hand, C2 = 0, we have at (xq, Aq) a "fold of order m" if A'^'^^(e) = for A; < m [12]. 
In the case of a simple singular point, we have either 

(i) dim{Af{Fu{uo, Xo))) = 1, and 
(ii) FA(Mo,Ao)G7^(F„(Mo,Ao)), 

or we have 

(i) dim(Ar(F,(no,Ao)))=2, and 
(ii) FA(Mo,Ao)^7^(F„(Mo,Ao)). 

In this case, we have a situation of branch-switching, in which there are are two 
branches of solutions crossing the point {uq, Ao) [12]. With high-dimensional null- 
spaces, the situation becomes increasingly complicated, with multiple branching 
solutions of various forms. We illustrate the three representative cases as they 
would appear in (HmH, A) space in Figure 1. 

We conclude, for completeness, this section with an exposition of the generalized 
form of the bifurcation analysis method of Lyapunov-Schmidt, as used for this 
problem by Walsh [:' ]. The exposition follows Zeidler [17]. We assume that 

(i) Fuiuf), Ao) is a Fredholm operator of index k 
(ii) dim(Ar(F„(Mo,Ao))) =n 

Then we define projection operators P : X — )■ X and Q : X -^ X with P{X) = 
Af{Fu{uo,Xo)) and (/ - Q)(Y) = 7^(F„(Mo, Ao)). Now the equation F{u,X) = is 
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Figure 1. Coinmon Locally Bifurcating Solution Paths. 



equivalent to the pair 



iI-Q)F{y + z,\) = 
QF{y + z,X) = 



(4.3) 
(4.4) 



with y = {I — P)u and z = Pu. Now the first equation satisfies the assumptions of 
the Implicit Function Theorem, and so we can get a unique solution branch y{z, A), 
substitute the solution in the second equation to obtain the branching equation 

QF{y{z,X) + zA) = ^ (4.5) 

then solve for z{\) to get the branch u = y{z{\), A) + z{\) Note that, in practice, 
one solves (4.3) by expanding the operators in the bases of J\f{Fu{uo, Xq)) and 
M{Fu{uo, Ao)*). Hence, to be constructive, one must already have an estimation of 
the dimension of the null-space of F„(mo, Aq). In particular, Walsh [ ] performs 
the analysis with the starting assumption that dim(A/'(F„(Mo, Aq))) = 1. So while 
the technique is appropriate once this is known, it is not a standalone method of 
bifurcation analysis. 

5. Numerical Bifurcation Theory 

In the procedure of continuation, we seek to find a solution m to a problem 
F{u, A) = as we move along A. In the case of a nonsingular Jacobian of the 
operator {Fu{uo, Aq)), it is standard to apply Newton's method on a discretization 
of the solution space. So we perform continuation by repeatedly iterating A = 
A + AA then resolving F[u, A) = for u. However, the procedure is invalid in 
the case of a singular Jacobian of the equation operator, necessitating alternatives 
for traversing a solution as the parameter varies. Depending on the form of the 
bifurcation, various procedures exist. Considering the finite-dimensional case (e.g. 
for a discretization of u), we have that F{u, A) maps M^ x M to M^. In this case, 
the rank of [Fu, Fx] = N ii either Fu is nonsingular, which is the case above, or 
Fx{uo, Ao) ^ 7l{Fu{uo, Aq)). In this case di'm{Af{[Fu, Fx])) = 1. Let's say we have a 
solution F{uo, Aq) = 0. In the latter case, we cannot just set Ai = Aq + A (with a 
constant A) and solve for ui, but we can solve for F{ui, Ai) = together with one 
more scalar equation, which we can set to be a constraint on the total magnitude in 
the change of (m. A). With pseudo-arclength continuation, at a certain parameter 
Ao and solution vector uq, and a direction vector [iiQ, Aq) of the solution branch 
determined thus far, you run Newton's method on the two equations F{ui, Ai) = 
and (mi — uo)*uo + (Ai — Ao)A — As = 0, where As is a constant "arclength" term. 
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It can be shown that the Newton's method Jacobian matrix is nonsingular if the 
point is either one at which F„ is nonsingular or a fold [6]. 

If diin{J\f{Fu)) > 1 or F\{uo, Aq) G 7l{Fu{uo, Aq)) then the Jacobian of Newton's 
method in the pseudo-arclength continuation equations is singular as well, and 
other procedures must be used to continue the solution. Methods of continuation 
usually involve constructive techniques, wherein the nullspace vectors for F^ and 
the solution F^cpr = -Pa are explicitly found and coefficients constructed (for a list 
of relevant algorithms, see [12]). With even higher-dimensional null-spaces, this 
framework is generalized. 

6. Setup for Hamiltonian Constraint Bifurcation 
We reprint the equations (2.2)-(2.3) here 

V^'ijj + 2Trpip^ = (6.1) 

^ = 0, r = (6.2) 

ip = l, r = l (6.3) 

and perform a reduction for bifurcation analysis following [6]. 

The goal is to transform this equation into a form wherein its linearization be- 
comes an eigenvalue problem. This can be done by parameterizing an extra bound- 
ary condition. The equation can be rewritten in an equivalent form as 

VV + 27rpV^ = (6.4) 

^{0)=p (6.5) 

^(0) = (6.6) 

Now we seek to solve F{p,p) = ^jj{l,p,p) — 1 = 0. Writing ipp{t,p,p) = -j^{t,p,p), 
Fp{p,p) = ipp{l,p,p) and similarly Fp{p, p) = ipp{l,p,p) Now ipp satisfies 

V^Vp + lOvrpV'Vp = (6.7) 

MO) = 1 (6-8) 



di!p 



(0) = (6.9) 



dr 

Which is an eigenvalue problem with a distinct solution. Define a{p, X) = 4'p{l, p, p) = 
Fp{p,p). Similarly V'p satisfies 

V^Vp + 27rV^^ + lOvrpV'Vp = (6.10) 

Mo) = o (6-11) 



dil^P 



(0) = (6.12) 



dr 
and likewise define b{p,p) = ipp(l,p,p) = Fp{p,p). Now we have -F'(p,p) = ia,b), So 

iV(F(pp))) can be I „ J , I , J , both, or neither. Let us discuss each of these cases 

separately. If the null-space is empty, then we have a nonsingular Jacobian of the 
operator, so continuation can proceed with Newton's method as standard. If the 

null-space is I -, ) , so b{p, p) = 0, then this is a situation where the same solution 
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exists for a differential increase in p, whicli is also an uninteresting case. If the null- 
space is I „ J , so a{p, p) = and b{p, p) ^ 0, we now have a case where N{Fp) = 1 

but Fp ^ Fp, so we have a fold. Pseudo-arclength continuation must be performed, 
and the form of the continuation gleamed from the coefficients. Finally, in the case 
of a full two-dimensional null-space, we have reached another solution branch. Here, 
there are two solution branches intersecting, and so we have a choice between three 
directions in the continuation of the solution. We now seek to perform this analysis 
numerically to locate any critical points for (6.1) and see one of these scenarios 
occurs. 



7. Numerical Results 

Recall that a continuation procedure attempts to solve F{u, A) for varying A by 
a series of successive iterations along discrete steps of A. The procedure begins at 
some Ao, finds the solution F{uo, Aq), finds the nullspace information of Fu{uo, Aq) 
and Fx{uo,Xo), then using this information performs a suitable continuation step 
to find the next solution F{ui,Xi). We used the continuation software AUTO to 
trace the solutions to (2.2) with boundary conditions (2.3). AUTO ([7]) performs 
numerical continuation on ODEs (notice that due to symmetry the differential 
equation becomes an ODE with just variable r). AUTO calculates the dimension 
of the null-space of the Jacobian of the differential operator and performs, as re- 
quired, Newton's method, pseudo-arclength continuation or explicit calculation of 
the null-space directions if it identifies higher-order bifurcations, as discussed in 
(5). It uses the bordering algorithm for the case of pseudo-arclength continuation 
(['"]) and otherwise more intricate linear algebraic algorithms for higher-dimensional 
bifurcations. 

Continuation reveals a quadratic fold at a value of p around p^ ~ 0.35. For 
p < Pc there are two solutions to (2.2), at p = pc, there is exactly one solution, and 
at p > pc there are no solutions. AUTO finds that at p = pc the nullspace of the 
discrete Jacobian differential operator has a dimension of one. In addition, it finds 
that the 2nd derivative information indicates that the continuation has a quadratic 
form. At all other points along the solution curve, the Jacobian is nonsingular. 
Now we let the exponent vary. Writing the equation as V'^ip + 2TTpip°' = 0, we 
investigate the continuation for different values of a. Knowing that the Laplacian 
operator on its own has an invertible Jacobian, we expect a fold to appear for some 
value of a. We find that indeed, with a > 1 a fold appears, the curvature of the 
solution continuation becoming sharper with increasing a. In Figure 3 we see the 
bifurcation diagram for a = 1, 1.25, 5, 10. We'd like to note that, in contrast to 
the previous numerical evidence for non-uniqueness [ ], rather than looking at the 
form of the solution curve and making an educated guess at the presence of a second 
solution curve, we confirm that there is one and only one additional solution curve 
which connects to the primary solution set in a quadratic curve and there are no 
additional branches or any solutions past pc- 

Now we look at two representative solutions from the two branches. Taking the 
two solutions at p = 0.2 for the original problem a = 5, we see We note that this 
confirms the results of [■">], as the solutions have the form of a Sobolev function and 
furthermore the conformal factor is considerably larger, suggesting a higher ADM 
energy, for the upper branch solutions. 
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Figure 2. Solution curve for (6.1). 



8. Discussion 

The Conformal Thin Sandwich approach has been a promising formulation of 
the Einstein constraints. The XCTS method appears to have the benefit of more 
physically intuitive and less computationally onerous specification of free data In 
contrast to its namesake predecessors, the conformal method and thin sandwich 
approach. In particular, you save the task of having to specify a divergence-free 
part of a symmetric tracefree tensor in the conformal method. As such it has been 
used extensively in the construction of the geometry of binary neutron stars and 
black holes (see the bibliography 5-12 in [ ]). As such this paper provides a point 
of caution for numerical relativists. 

As Baumgarte et al. [>] noted, the upper branch corresponds to a higher level 
of ADM energy than the lower branch in this model of a constant density black 
hole. As such, in most applications with this set of asymptotically fiat, spherically 
symmetric free data, the lower branch should be chosen as the correct, physically 
representative solution, but numerical relativists should be aware that a solver could 
reach either one and be able to identify the need to switch to the other branch. 
In addition, as p —)■ pc standard numerical procedures will become increasingly 
ill-conditioned. Pinning the solution to the expected branch is recommended as a 
solution strategy. 

Finally, it should be noted that since this bifurcation analysis is complete, rather 
than simply constructive, and so demonstrates conclusively that 1) there are exactly 
two solutions for p < pc and 2) there are no solutions for p > pc for this choice of 
free data. The fact that the same choice of free data is consistent with two qualita- 
tively disparate geometries suggests future investigation in the relation between the 
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Figure 3. Bifurcation diagram for a = 1, 1.25, 5, 10. 

structure in the free and dependent set of data. In addition, it would be an inter- 
esting investigation as to why a higher density is inconsistent with time-symmetry, 
conformal flatness and spherical symmetry. 



9. Conclusion 

In this article, we have examined the apparent bifurcation phenomena in the 
XTCS formulation of the Einstein constraints in a methodical way, using modern 
techniques in bifurcation theory and in numerical homotopy methods. We first 
gave an overview of the Einstein constraints in §2, and followed this in §3 with 
a brief introduction to the mathematical foundations for, and in §4 the frame- 
work for analyzing bifurcation phenomena in nonlinear operators equations. In §5, 
we developed the main ideas behind the construction of numerical homotopy, or 
path-following, methods in the numerical treatment of bifurcation phenomena, and 
in §6 we presented the set up of and in §7 we applied the continuation software 
package AUTO to the constraint problem. We verified the presence of the fold 
with homotopy-based numerical methods, confirming the earlier results. Analyz- 
ing the Hamiltonian constraint for time-symmetric conformally fiat initial data, as 
in [16, 5], we demonstrated the existence and location of a critical point, evidence 
that the solution branch at the critical point is one-dimensional, and the form of 
the solution as continued past the critical point is a simple quadratic fold. We con- 
firm Walsh's [lU] constructive Lyapunov-Schmidt analysis by showing numerically 
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Figure 4. Solutions for a = 5, p = 0.2. 



that there is indeed a solution point at which there is a one- dimensional kernel 
for the linearization of the differential equation operator, justifying the assumption 
made in his analysis, as well as confirming numerically that indeed the form of the 
solution continuation expansion coefficients are such that this equation exhibits a 
quadratic fold rather than branching or higher-order folds. 

The techniques presented here can be viewed as providing a framework for a 
more careful exploration of the solution theory for various parameterizations of the 
constraint equations, as well as the geometric relationship between the free data in 
the model problem investigated. 
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